Expanding or shrinking? range shifts in wild ungulates under climate change in Pamir-Karakoram mountains, Pakistan

Climate change is expected to impact a large number of organisms in many ecosystems, including several threatened mammals. A better understanding of climate impacts on species can make conservation efforts more effective. The Himalayan ibex (Capra ibex sibirica) and blue sheep (Pseudois nayaur) are economically important wild ungulates in northern Pakistan because they are sought-after hunting trophies. However, both species are threatened due to several human-induced factors, and these factors are expected to aggravate under changing climate in the High Himalayas. In this study, we investigated populations of ibex and blue sheep in the Pamir-Karakoram mountains in order to (i) update and validate their geographical distributions through empirical data; (ii) understand range shifts under climate change scenarios; and (iii) predict future habitats to aid long-term conservation planning. Presence records of target species were collected through camera trapping and sightings in the field. We constructed Maximum Entropy (MaxEnt) model on presence record and six key climatic variables to predict the current and future distributions of ibex and blue sheep. Two representative concentration pathways (4.5 and 8.5) and two-time projections (2050 and 2070) were used for future range predictions. Our results indicated that ca. 37% and 9% of the total study area (Gilgit-Baltistan) was suitable under current climatic conditions for Himalayan ibex and blue sheep, respectively. Annual mean precipitation was a key determinant of suitable habitat for both ungulate species. Under changing climate scenarios, both species will lose a significant part of their habitats, particularly in the Himalayan and Hindu Kush ranges. The Pamir-Karakoram ranges will serve as climate refugia for both species. This area shall remain focus of future conservation efforts to protect Pakistan’s mountain ungulates.

Introduction (SDMs) which are widely adopted in investigations of species distribution and range shifts [36,37].

Collection of presence records
Himalayan ibex and blue sheep presence records (Fig 1) were collected using two methods: camera trapping and double observer surveys. Cameras were left operational for 10 days in the first camera trapping in KNP, but in the latter surveys they were left operational for 40 days to increase the capture rate [42,43].
2. Double observer Survey: We carried out this survey in 2012-2016 in different parts (KNP, Gojal Valley, Shigar Valley, in Skardu district, and in Gilgit district) of the study area by dividing it into smaller blocks based on watersheds. These watersheds were not larger than daily ungulate/human movement ability. Two observers were sent for survey separated by time (15 minutes) if only one trail was available, or by space, if two trails were available. Each watershed was surveyed by walking along pre-determined routes [44]. The locations where Himalayan ibex and blue sheep were sighted, have been used as presence points to build the MaxEnt model.
We collected 143 and 60 presence points for Himalayan ibex and blue sheep, respectively (S1A and S1B Fig). We then screened these presence points in ArcGIS 10.7 (ESRI, Redland, USA) using nearest neighbor analysis to check spatial autocorrelation [36,43,45]. This analysis revealed a high clustering among presence points. Aggregation was, therefore, spatially filtered using SDMTools [46] to ensure independence [36,43,47]. This operation led to 36 and 29 presence points for Himalayan ibex and blue sheep, respectively, which we used in MaxEnt models (Fig 1).

Climatic variables
We downloaded 19 climatic variables from WorldClim 1.4 (https://www.worldclim.org/ current) [48] to predict currently suitable areas for Himalayan ibex and blue sheep. All the variables were in raster files (grid) with 30-arc second resolution (0.93 × 0.93 km = 0.86 km 2 at the equator). Further details and information on the realization and interpretation of the WorldClim variables used in this study can be found at https://pubs.usgs.gov/ds/691/. We checked all variables for multicollinearity and excluded highly correlated variables i.e., r � 0.70 (Pearson's correlation coefficient) [43]. This process led to use in the modeling analysis of six environmental variables: annual mean temperature (C˚), mean diurnal range (˚C), temperature seasonality [(standard deviation � 100) (˚C)], mean temperature of wettest quarter (˚C), mean precipitation (mm), and precipitation seasonality (%).
We used global circulation models (GCMs) MIROC5, BCC-CSM1-1, CCSM4, and Had-GEM2ES to predict the future distribution of Himalayan ibex and blue sheep under climate change conditions. Various organizations developed these models under the Coupled Model Intercomparison Project, phase 5 (CMIP5) and are considered highly reliable [36,49]. The future projections of these GCMs are based on representative concentration pathways (RCPs) which are greenhouse gas (GHG) concentration trajectories on a range of radiative forces suggested in the Intergovernmental Panel on Climate Change's (IPCC) fifth assessment report [50]. We used RCP 4.5 and RCP 8.5 the former is a moderate GHG mitigation scenario [51]  where emissions will peak around 2040 and then decline, while the latter is a scenario where GHG emissions will be the highest of all four RCPs (2.6, 4.5. 6.0 and 8.5) throughout the 21st century [27].

Modeling procedure
We used MaxEnt ver. 3.4.1 [52] to predict the current and future distribution of C. ibex sibirica and P. nayaur in Pakistan [25]. MaxEnt is a piece of machine learning software used to develop SDMs [53][54][55]. It is capable of predicting species distribution using presence-only data [56] and predicting the distribution of poorly known species [36,57]. We built the model using a logistic output format to yield environmental suitability ranging from 0 (unsuitable) to 1 (highly suitable) [58]. We fixed the regularization multiplier to 1, selected 5,000 iterations [27], and ran 20 replicates with cross-validations tests [43].
Different GCM projections can have inherited uncertainties. To avoid this, we used area under the curve (AUC) scores as weighting coefficients that resulted from 20 cross-validations for each of four GCMs and produced a single forecast for each time scale by averaging all individual GCMs for that time slice. [28,[59][60][61]. We used ten percentile training presence values as the threshold to develop binary presence/absence maps [43].
The model was projected to entire GB. To project the models calibrated for survey area over entire GB, the variables in the projection area must meet a condition of environmental similarity with the environmental data used for calibrating the model. Therefore, we preliminarily ascertained that this condition was verified for both current and future projections by inspecting Multivariate Environmental Similarity Surfaces (MESS), the MESS calculates the similarity of each point in the region of projection to a set of reference points (e.g., background data) and maps the results [56] MESS maps produced by MaxEnt can help users identify extrapolated areas and provide a quantitative measure of projection uncertainty.

Model validation
We tested the predictive performance of the models with different methods: receiver operated characteristics, analyzing the AUC [62], and the true skill statistic (TSS) [63]. AUC assesses models' discrimination ability with values ranging from 0 (equaling random distribution) to 1 (perfect prediction). TSS compares the number of correct forecasts minus those attributable to random guessing, to that of a hypothetical set of perfect forecasts. It considers both omission and commission errors and success as a result of random guessing. Its values range from -1 (a performance no better than random) to +1 (perfect agreement).

Niche overlap
We calculated the niche overlap between C. ibex sibirica and P. nayaur for predicted habitats using ENMTools [64] in the current time and future climate change scenarios. ENMTools uses MaxEnt map values of habitat suitability for each grid and measures niche overlap using D and I values [64]. It uses Schoener's D value to calculate niche overlap and gives probability distributions with values ranging from 0 (no overlap) to 1 (complete overlap). Similarly, Hellinger's I-statistic in ENMTools measures models' ability to estimate true suitability [64].

Model performance
The AUC values for our models were 0.969 ± 0.025 and 0.821 ± 0.138 for blue sheep and Himalayan ibex, respectively. TSS values were 0.841 ± 0.007 and 0.454 ± 0.281 for blue sheep and Himalayan ibex, respectively. Both tests suggest strong performances of our models.

Current distribution of Himalayan ibex and blue sheep
Our binary maps showed ca. 26 500 km 2 (37.71% of total study area) and ca. 6 500 km 2 (9.26% of total study area) suitable for Himalayan ibex and blue sheep, respectively (Fig 2).
We found that the current habitat predicted for Himalayan ibex included the latitudes from 34˚to 37˚and the longitudes from 73˚to 77˚. The most suitable habitats fell in the Karakoram range, followed by the Hindu Kush, and then to a minor extent in the Himalayas (Fig 2A). The habitat suitability of Himalayan ibex was predicted in all ten districts of GB with strongholds in Hunza, Nagar, Shigar, and Ghanche districts. We found that habitats suitable to blue sheep were between the latitudes 35˚to 37˚and the longitudes 74˚to 77˚along the Pakistan-China border in the Pamir-Karakorum range that administratively falls in Hunza district, followed by some parts of the Shigar and Ghanche districts along the Pakistan-China border (Fig 2B). We found that annual mean precipitation, mean temperature of the wettest quarter, and temperature seasonality were the most important variables (with 91.6% contribution) in predicting suitable habitats for blue sheep (S1 Table). For ibex, annual mean precipitation, annual mean temperature, and precipitation seasonality were key habitat predictors with an 89% contribution (S2 Table). In the extreme climate change scenario (RCP 8.5 of 2070), blue sheep lost (58%) from the suitable areas that it has currently occupied and gained new suitable areas by extending its current range towards the east. Himalayan ibex gained the least and lost (64.80%) in RCP 8.5 of 2070 (Table 3 and Figs 5 and 6). The model predicted habitat shrinkage to an area of 2,515 km 2 for blue sheep and 9,248 km 2 for ibex under the extreme climate change scenario.
The center of suitable Himalayan ibex habitat gradually shifted from the north to the east in RCP 4.5 and RCP 8.5 of 2050, and RCP 4.5 of 2070, while in RCP 8.5 of 2070, it again shifted from the east to the north. The center of the suitable habitat of blue sheep first shifted gradually from the west towards the north in RCP 4.5 and RCP 8.5 of 2050, and RCP 4.5 of 2070. In RCP 8.5 of 2070, it shifted towards the east from the north. The MESS analysis predicted some areas with novel climate conditions across the range for both P. nayaur and C. ibex sibirica in the future projections. However, these areas were found outside the training range of our model (S1-S8 Figs).  (Table 4 and Fig 7).

Discussion
The use of SDMs for the predictive distribution of biodiversity [65] has increased as the approach is considered efficient in predicting species distribution and climate change impact [66] which aids in species conservation planning [55]. MaxEnt is widely used for its proven ability to construct models using presence-only data [67]. This model worked well on our presence data as indicated by the AUC scores (>0.8), which places it among the best-published models [25,26,28,68]. The higher TSS values further supported the credibility of results [36,47]. The melting of Himalayan glaciers has increased in the 21 st century [69] while the glaciers of the Hindu Kush and Karakoram will melt at a slower rate [70]. In fact, some glaciers in the  higher watersheds of the Karakoram are expanding [71] although at the same time they are thinning. However, regardless of the three described scenarios, the snow on these glaciers regulates ecological processes and patterns [72] and any change in glacier mass, negative or positive, will affect associated biodiversity. Our results for habitat loss and gain were strikingly aligned with the existing knowledge on glaciology. We found that global climate change will have significant effects on the habitats of mountain ungulates in northern Pakistan, though these effects are more pronounced in Hindu Kush, and Himalaya ranges. Our model for current time predicted 6,510 km 2 and 26,510 km 2 of suitable area for blue sheep and Himalayan ibex, respectively. Both model species are present in most of the predicted habitats, or they occupied those areas historically [30,33]. Ironically, Khan et al., (2014) reported sighting records of ibex in Tangir Valley of Diamer district, which is beyond the suitable habitat predicted in the current study, as well as outside of the former IUCN range [73]. This probably indicates southwards expansion of ibex in recent years. Our model predicted suitable habitat for blue sheep on the Braldu glacier where sheep do not currently exist [74]. Interestingly, older records indicate the presence of blue sheep in this area, e.g., [29] quote a sighting by T. J. Roberts in this area in 1975.
Both blue sheep and Himalayan ibex habitats are usually between the timber and snow lines at elevations of 3,500-5,500 m, and differ as Blue sheep prefers habitats with steep rolling hills and Himalayan ibex prefer precipitous habitats [33]. These habitats are usually devoid of thick vegetation. Hence, precipitation is a vital factor to sustain life in this zone. We found annual precipitation to be the most contributing variable in predicting suitable habitat for both blue sheep and Himalayan ibex. Annual mean temperature was the second most important variable for Himalayan ibex, and temperature of wettest quarter the second most important for blue sheep. The dry habitats of both ibex and blue sheep have short growing seasons, and any weather fluctuation might leave species starving [75]. Artemisia and Ephedra shrubs are described as the ibex's main food sources [33]. A year of good winter precipitation and normal mean summer temperature enables shrubs to maximize their growth and green cover [76]. Blue sheep's preferred diet comprises of grasses, forbs, and shrubs Berberis, Polygonum, and Ephedra, respectively [33]. Even in the summers, precipitation at elevations above 4,000 m can bring temperatures below zero and constraint vegetative growth [76]. Hence, temperatures of wettest quarters (June, July, and August) play a decisive role in selecting suitable habitat for blue sheep. Khan et al. (2016) found annual precipitation and minimum temperature to be important variables for developing suitability models for C. ibex sibirica and P. nayaur, respectively. Aryal et al. (2016) and Luo et al. (2015) reported annual mean temperature as the most influencing variable in predicting suitable habitat for P. nayaur.
We observed a sharp decline (56% in RCP 4.5 and 58% in RCP 8.5) in the currently available suitable habitat for blue sheep and (33.70% in RCP 4.5 and 64.80% in RCP 8.5) for Himalayan ibex in extreme climate change scenarios for 2070. This is consistent with [25]who observed a decrease in blue sheep suitable habitat in the future due to climate change in Nepal. Similarly, Luo et al. (2015) reported a 30-50% range reduction for ungulates on the Tibetan plateau under different climate change scenarios.
Climate drives evolutionary processes, forcing animals to migrate to higher elevations or extend their distributional ranges towards the Northern Hemisphere [77] or eastward direction [28]. This process is believed to have occurred in the Miocene Epoch when members of the Caprinae in Eurasia and Africa began inhabiting the newly formed mountain ranges of the Himalayas, Karakoram, Hindu Kush, and Pamirs, which emerged from the sea during the Tertiary Period [33]. We expect a similar migration in northern Pakistan because the centers of predicted suitable habitat for Himalayan ibex will shift from north to east in RCP 4.5 and RCP 8.5 of 2050 and 2070 and again from east to the north in RCP 8.5 of 2070. For Himalayan ibex, it will shift from west to north in RCP 4.5 and RCP 8.5 of 2050 and 2070 and from north to east in RCP 8.5 of 2070. Species co-evolved over millions of years, enabling them to co-exist by selecting different niches [78]. Our model predicted a moderate niche overlap between blue sheep and Himalayan ibex, and this overlap was predicted to increase if the extreme climatic conditions assumed in future scenarios prevail. Increasing temperatures and precipitation have already impacted Himalayan flora [79]. Alpine habitats have short growing seasons [80,81] and offer relatively few species of grasses, sedges, forbs, shrubs, ferns, lichens, and mosses to Himalayan ibex and blue sheep [82][83][84]. Hence, these climatic changes in alpine ranges will increase the chances of habitat mismatch for many floral species [28,80]. Climate change, together with anthropogenic effects transforming land for agriculture or afforestation, road construction, and mining could further shrink habitats suitable for ungulates [28,68], potentially affecting their perpetuity and the proper functioning of ecosystems [85,86].
Conservationists emphasize on locating habitats likely to be least affected by climate change and continue serving as suitable habitats (future refugia), and protecting them from anthropogenic activities [21,87,88]. Our model predicted such climate refugia for Himalayan ibex to be comprised of three national parks: Khunjerab National Park (KNP), Central Karakoram National Park (CKNP), and Qurumbar National Park (QNP) (Fig 6). For blue sheep, such refugia exists in the buffer zone of KNP, along with a few patches on the Braldu glacier of CKNP (Fig 5). It is noteworthy, however, that Himalayan ibex will lose most of its current suitable habitat in CKNP in Baltistan division and areas around QNP in the future, but the areas of CKNP in Nagar district will remain stable. All three mountain ranges in our study area provide vital habitats to several mountain ungulates. Unfortunately, most of suitable habitats in Hindu Kush and Himalayas are expected to be altered under future scenarios. On contrary, the Pamir-Karakoram is likely to remain stable and continue accommodating both Himalayan ibex and blue sheep. The relatively lower effect of climate change in this range is likely due to the barrier effect of the Hindu Kush and Himalayas which blunt the monsoon, helping maintain the aridity of the Karakorum's' alpine steppes [21,71].

Conclusions
Our study demonstrate that the current suitable habitat of Himalayan ibex and blue sheep are vulnerable to climate change. Under the rapid climate change Himalayan ibex will lose most of its current suitable habitat in Himalayans and Hindu Kush while blue sheep that currently exists only in Pamir-Karakoram range will be slightly affected. The current network of protected areas (KNP and CKNP) will serve climate refugia for mountain ungulates.
There is urgent need to revisit protected areas management strategies in Pakistan, to enhance their effectiveness for conservation of mountain ungulates. The finding of this study can be used to revisit or align boundaries of existing protected areas with the future predicted habitats. Management and protection efforts shall remain disproportionally higher in parks that encompass climate refugia for mountain ungulates of the region.